model { for(i in 1:N){ mu_M[i] <- alpha*x[i] M[i] ~ dnorm(mu_M[i],prec1) mu_y[i]<- c*x[i]+beta*M[i] y[i] ~ dnorm(mu_y[i],prec2) j1[i]~dunif(0.5,100.5) j2[i]~dunif(0.5,100.5) j3[i]<- round(j1[i]) j4[i]<- round(j2[i]) M1[i] ~ dnorm(alpha,prec1) M0[i] ~ dnorm(0,prec1) mu_y1[i]<- c+beta*M1[i] mu_y0[i]<- beta*M0[i] te[i]<-(mu_y1[i]-mu_y0[i])/deltax mu_M2[i] <- alpha*(x[j3[i]]) M2[i] ~ dnorm(mu_M2[i],prec1) mu_y2[i]<- beta*M2[i] mu_M3[i] <- alpha*(x[j4[i]]) M3[i] ~ dnorm(mu_M3[i],prec1) mu_y3[i]<- c+beta*M3[i] de[i]<-(mu_y3[i]-mu_y2[i])/deltax ie[i]<-te[i]-de[i] } alpha ~ dnorm(0.0,0.000001) beta ~ dnorm(0.0,0.000001) c ~ dnorm(0.0,0.000001) var1 ~ dgamma(1,0.1) var2 ~ dgamma(1,0.1) prec1 <-1/var1 prec2 <-1/var2 }